%% it produces Figure 13


get_params;
% load calibration South at initial steady state
str_region="South_";
str1="params_";
str_3="calibration_Table2";
str_save=append(path_calibration,str1,str_region,str_3,".mat");
filename=cell2mat(str_save);
load(filename)
read_pars;read_pars_s;

k_vec   = 1.0:-0.025:0.3;
warning 'off'
psi   =-log(.5)/10;
dt    = 1;
T_y   = 100;
T_sim = T_y/dt;

inde_samechi=0;


%% load calibration South at initial steady state
str_region="South_";
str1="params_";
str_3="calibration_Table2";
str_save=append(path_calibration,str1,str_region,str_3,".mat");
filename=cell2mat(str_save);
load(filename)
read_pars;read_pars_s;
Y=1;
delta_s = .0888;%exit prob
Profit  =A_s/nu;
cost_s=chi_s*A_s;
mu_s    = Profit-chi_s*A_s;
crit_I=1;
while crit_I>1/100000000
    v0_s    = fzero(@(v) (r+delta_s)*v - theta/(1+theta)*zeta*v.^(1+1/theta)-mu_s,10);
    z0_s    = z_0_s;
    Avg_z_s   = zavg;
    iota_s    = delta_s*(1-z0_s/Avg_z_s);
    zeta_s= iota_s* v0_s^(-1/theta);

    k_s_ben=  v0_s*z0_s;
    m_s_ss  = (Y/A_s)^(nu/(nu-1))/(Y*Avg_z_s);
    crit_I =abs(zeta-zeta_s);
    zeta=zeta_s;
end


% load calibration North at initial steady state
str_region="North_";
str1="params_";
str_3="calibration_Table2";
str_save=append(path_calibration,str1,str_region,str_3,".mat");
filename=cell2mat(str_save);
load(filename)
read_pars;read_pars_n;



k_n_ben=k;

sigma_n =sigma;
delta_n=0.0775;
Profit=A_n/nu;
cost_n=chi_n*A_n;
mu_n    = Profit-chi_n*A_n;
crit_I=1;
while crit_I>1/1000000
    [v0_n,fv]    = fzero(@(v) (r+delta_n)*v - theta/(1+theta)*zeta*v.^(1+1/theta)-mu_n,1);
    z0_n     = z_0_n;
    Avg_z_n  = zavg;
    iota_n    = delta_n*(1-z0_n/Avg_z_n);
    zeta_n= iota_n* v0_n^(-1/theta);
    k_n_ben  =  v0_n*z0_n;
    m_n_ss    = (Y/A_n)^(nu/(nu-1))/(Y*Avg_z_n);
    crit_I    = abs(zeta-zeta_n);
    zeta  = .85*zeta_n+.15*zeta;
end


l_s_0=1;
l_n_0=1;
dt   = 1;
disc = (1/(1+r)).^(0:dt:T_y-dt);
l_reallocation = 1;
demand_ext     =1;
inve_ext=1;


W_vec=[];
W_vec_s=[];
W_vec_n=[];

W_vec=[];
W_vec_s=[];
W_vec_n=[];
for ip=1:length(k_vec)
    k_n = k_n_ben;
    k_s = k_vec(ip)*k_s_ben;

    v0_s      = k_s/z_0_s;%new value of equity
    mu_s_new  = v0_s*(r+delta_s) - theta/(1+theta)*zeta_s*v0_s^(1+1/theta);
    A_s_new   = (mu_s_new+cost_s)*nu;
    iota_s_new= zeta_s* v0_s^(1/theta);
    if iota_s_new>delta_s,error('exploding path'),end

    v0_n      = k_n/z_0_n;
    mu_n_new  = v0_n*(r+delta_n) - theta/(1+theta)*zeta_n*v0_n^(1+1/theta);
    A_n_new   = (mu_n_new+cost_n)*nu;
    iota_n_new= zeta_n*v0_n^(1/theta);
    if iota_n_new>delta_n,error('exploding path'),end

    %% initial steady state
    Y_ss=1;l_n_ss=1;l_s_ss=1;l_s=1;l_n=1;
    m_s_0       = m_s_ss;
    m_n_0       = m_n_ss;

    %% Find new steady state
    Y_new = 1;crit_Y=1;
    if demand_ext==1
        while crit_Y>1/100000
            w_s       = (nu-1)/nu*(A_s_new/Y_new)^(1/(1-nu));
            w_n       = (nu-1)/nu*(A_n_new/Y_new)^(1/(1-nu));
            if l_reallocation==1
                [l_n,fv]      =  fsolve(@(ell) w_n - w_s - h_n*max(0,ell).^eta  + h_s*max(0,((1-ell*wei_n)/wei_s)).^eta,1,options);
                if abs(fv)>1/10000,error('fv'),end
                %if l_n>1 || l_n<0, error('l_n'),end
                l_s           = (1-l_n*wei_n)/wei_s;%l_s*wei_s+l_n*wei_n=1
            end
            Y_imp         =  (wei_s*l_s*A_s_new^(1/(1-nu)) + wei_n*l_n*A_n_new^(1/(1-nu)))^((nu-1)/(nu-2));
            crit_Y        = abs(Y_imp- Y_new) ;
            Y_new         = Y_new*(1+0.1*(Y_imp- Y_new)/Y_new);
        end
    else
        w_s       = (nu-1)/nu*(A_s_new/Y_new)^(1/(1-nu));
        w_n       = (nu-1)/nu*(A_n_new/Y_new)^(1/(1-nu));
        if l_reallocation==1
            [l_n,fv]  =  fsolve(@(ell) w_n - w_s - h_n*max(0,ell).^eta  + h_s*max(0,((1-ell*wei_n)/wei_s)).^eta,1,options);
            l_s       = (1-l_n*wei_n)/wei_s;%l_s*wei_s+l_n*wei_n=1
        end
    end
    l_n_ss_new     = l_n;
    l_s_ss_new     = l_s;

    %% Equilibrium across regions: output
    time=0;inde_sim=0;l_n=1;l_s=1;Y_new = 1;
    z_tm1_s=Avg_z_s*m_s_ss;mm1_s=m_s_ss;  z_tm1_n=Avg_z_n*m_n_ss;mm1_n=m_n_ss;
    while time<T_y
        time=time+dt;inde_sim=inde_sim+1;
        if demand_ext==1
            Y_new         =  (wei_s*l_s*A_s_new^(1/(1-nu)) + wei_n*l_n*A_n_new^(1/(1-nu)))^((nu-1)/(nu-2));
        else
            Y_new         = 1;
        end

        w_s       = (nu-1)/nu*(A_s_new/Y_new)^(1/(1-nu));
        w_n       = (nu-1)/nu*(A_n_new/Y_new)^(1/(1-nu));

        %find the average productivity in the transition
        m0_new_s=delta_s*m_s_ss;m0_new_n=delta_n*m_n_ss;
        crit_m=1;
        while crit_m>1/1000000
            m_s_0    = mm1_s*(1-delta_s)+m0_new_s;
            m_n_0    = mm1_n*(1-delta_n)+m0_new_n;
            z_t_s    = z_tm1_s +z_tm1_s*(iota_s_new-delta_s) +  z0_s*m0_new_s;
            z_t_n    = z_tm1_n +z_tm1_n*(iota_n_new-delta_n) +  z0_n*m0_new_n;
            z0_t_s=z_t_s/m_s_0;
            z0_t_n=z_t_n/m_n_0;
            m_s       = (Y_new/A_s_new)^(nu/(nu-1))/(Y_new*z0_t_s);
            m_n       = (Y_new/A_n_new)^(nu/(nu-1))/(Y_new*z0_t_n);
            m0_new_s_impl=m_s-mm1_s*(1-delta_s);
            m0_new_n_impl=m_n-mm1_n*(1-delta_n);
            crit_m=max(abs(m0_new_s_impl-m0_new_s),abs(m0_new_n_impl-m0_new_n));
            m0_new_s=m0_new_s_impl; m0_new_n=m0_new_n_impl;
        end
        %update the past values
        z_tm1_s=z_t_s;z_tm1_n=z_t_n;mm1_s=m_s;mm1_n=m_n;

        y_s       = A_s_new*m_s*z0_t_s;
        y_n       = A_n_new*m_n*z0_t_n;
        if inde_sim==1
            i_s       = (m_s_0*delta_s + m_s-m_s_0)*k_s_ben + 1/(1+theta)*m_s*z0_t_s*zeta_s*(iota_s_new/zeta_s)^(1+theta);
            i_n       = (m_n_0*delta_n + m_n-m_n_0)*k_n_ben + 1/(1+theta)*m_n*z0_t_n*zeta_n*(iota_n_new/zeta_n)^(1+theta);
        else
            i_s       = (m_s*delta_s)*k_s_ben+ m_s*iota_s_new*z0_t_s;
            i_n       = (m_n*delta_n)*k_n_ben+ m_n*iota_n_new*z0_t_n;
        end
        W_time_s_vec_2(inde_sim,ip) = (y_s-m_s*cost_s*z0_t_s -i_s -h_s*l_s^(eta));
        W_time_n_vec_2(inde_sim,ip) = (y_n-m_n*cost_n*z0_t_n -i_n -h_n*l_n^(eta));
        W_time_vec_2(inde_sim,ip)   = W_time_s_vec_2(inde_sim,ip)*l_s*wei_s+W_time_n_vec_2(inde_sim,ip)*l_n*wei_n;
        m_s_0=m_s;m_n_0=m_n;

        if ip==1 && inde_sim==1
            c_ss_dixit = (y_s-i_s)*l_s*wei_s+(y_n-i_n)*l_n*wei_n;
        end

        if  l_reallocation==1 &&  demand_ext==1
            cost_firm_s=m_s*cost_s*z0_t_s;
            cost_firm_n=m_n*cost_n*z0_t_n;
            W_time_s_vec(inde_sim,ip)     = y_s-i_s-cost_firm_s-h_s*l_s^(eta)*dt;
            W_time_n_vec(inde_sim,ip)     = y_n-i_n-cost_firm_n-h_n*l_n^(eta)*dt;
            W_time_vec(inde_sim,ip)       = W_time_s_vec(inde_sim,ip)*l_s*wei_s+W_time_n_vec(inde_sim,ip)*l_n*wei_n;
            y_time_mat(inde_sim,ip)       = y_s*l_s*wei_s+y_n*l_n*wei_n;
            y_time_s_mat(inde_sim,ip)     = y_s;
            y_time_n_mat(inde_sim,ip)     = y_n;
            i_time_mat(inde_sim,ip)       = i_s*l_s*wei_s+i_n*l_n*wei_n;
            i_time_s_mat(inde_sim,ip)     = i_s;
            i_time_n_mat(inde_sim,ip)     = i_n;
            cost_time_mat(inde_sim,ip)    = cost_firm_s*l_s*wei_s+cost_firm_n*l_n*wei_n;
            cost_time_s_mat(inde_sim,ip)  = cost_firm_s;
            cost_time_n_mat(inde_sim,ip)  = cost_firm_n;
            conj_time_mat(inde_sim,ip)    = h_n*l_n^(eta)*dt*l_n*wei_n+h_s*l_s^(eta)*dt*l_s*wei_s;
            conj_time_n_mat(inde_sim,ip)  = h_n*l_n^(eta)*dt;
            conj_time_s_mat(inde_sim,ip)  = h_s*l_s^(eta)*dt;
            l_time_s_mat(inde_sim,ip)     = l_s;
            l_time_n_mat(inde_sim,ip)     = l_n;
        end

        if l_reallocation==1
            u_s      = w_s - h_s*l_s^(1+eta);
            u_n      = w_n - h_n*l_n^(1+eta);
            dl_s     = l_s_ss_new-l_s_ss;
            l_s_new  = 1 + dl_s*(1-exp(-psi*time));
            l_n_new  = (1-l_s_new*wei_s)/wei_n;
            l_n      = l_n_new;
            l_s      = l_s_new;
        else
            l_n     = 1;
            l_s     = 1;
        end

    end

    welf_disc_all(ip)  =r/(1+r)*((disc*W_time_vec_2(:,ip))+disc(end)/r*W_time_vec_2(end,ip))

    if   l_reallocation==1 &&  demand_ext==1
        welf_disc(ip)      =r/(1+r)*((disc*W_time_vec(:,ip))+disc(end)/r*W_time_vec(end,ip));
        welf_disc_s(ip)    =r/(1+r)*((disc*W_time_s_vec(:,ip))+disc(end)/r*W_time_s_vec(end,ip));
        welf_disc_n(ip)    =r/(1+r)*((disc*W_time_n_vec(:,ip))+disc(end)/r*W_time_n_vec(end,ip));
        welf_disc_ss(ip)   =W_time_vec(end,ip);
        welf_disc_s_ss(ip) =W_time_s_vec(end,ip);
        welf_disc_n_ss(ip) =W_time_n_vec(end,ip);
        y_disc(ip)         =r/(1+r)*((disc*y_time_mat(:,ip))+disc(end)/r*y_time_mat(end,ip));
        y_disc_s(ip)       =r/(1+r)*((disc*y_time_s_mat(:,ip))+disc(end)/r*y_time_s_mat(end,ip));
        y_disc_n(ip)       =r/(1+r)*((disc*y_time_n_mat(:,ip))+disc(end)/r*y_time_n_mat(end,ip));
        i_disc(ip)         =r/(1+r)*((disc*i_time_mat(:,ip))+disc(end)/r*i_time_mat(end,ip));
        i_disc_s(ip)       =r/(1+r)*((disc*i_time_s_mat(:,ip))+disc(end)/r*i_time_s_mat(end,ip));
        i_disc_n(ip)       =r/(1+r)*((disc*i_time_n_mat(:,ip))+disc(end)/r*i_time_n_mat(end,ip));
        cost_disc(ip)      =r/(1+r)*((disc*cost_time_mat(:,ip))+disc(end)/r*cost_time_mat(end,ip));
        cost_disc_s(ip)    =r/(1+r)*((disc*cost_time_s_mat(:,ip))+disc(end)/r*cost_time_s_mat(end,ip));
        cost_disc_n(ip)    =r/(1+r)*((disc*cost_time_n_mat(:,ip))+disc(end)/r*cost_time_n_mat(end,ip));
        conj_disc(ip)      =r/(1+r)*((disc*conj_time_mat(:,ip))+disc(end)/r*conj_time_mat(end,ip));
        conj_disc_s(ip)    =r/(1+r)*((disc*conj_time_s_mat(:,ip))+disc(end)/r*conj_time_s_mat(end,ip));
        conj_disc_n(ip)    =r/(1+r)*((disc*conj_time_n_mat(:,ip))+disc(end)/r*conj_time_n_mat(end,ip));
    end
end
W_mat   =welf_disc_all';
%%
W_baseline_dixit= W_mat;
k_vec_dixit=k_vec;

[mm,index_max_coun]=max(W_baseline_dixit);
load subsidy_opt.mat
W_opt = subsidy_opt.W_opt;
cons_ss= subsidy_opt.cons_ss;
k_vec=subsidy_opt.k_vec_opt;
[mm,index_max_base]=max(W_opt);

figure(102)
plot(100*(1-k_vec_dixit),100/c_ss_dixit*(W_baseline_dixit)-100/c_ss_dixit*W_baseline_dixit(1),'--r','Linewidth',6),hold on
plot(100*(1-k_vec),100*(W_opt-W_opt(1)),'-b','Linewidth',6),hold off
xline(100*(1-k_vec(index_max_base)),':k','Optimal size: baseline','LineWidth',2,'FontSize',18,'Interpreter','Latex','LabelVerticalAlignment','bottom')
xline(100*(1-k_vec_dixit(index_max_coun)),':k','Optimal size: Dixit-Stiglitz','LineWidth',2,'FontSize',18,'Interpreter','Latex','LabelVerticalAlignment','bottom')
%axis([0 70 0 3])
le=legend('Dixit-Stiglitz','Baseline');
set(le,'Interpreter','Latex','Fontsize',16,'Location','northwest');
xlabel('$\vartheta+\tau$, \%','Fontsize',18,'Interpreter','Latex')
ylabel('Consumption gains, in \% dev. from  status quo','Fontsize',18,'Interpreter','Latex')
grid on
str1="/figures/figure_wefare_decomp_3.pdf";
str_save1=append(root,str1);
settt
print(gcf,'-dpdf',str_save1)


